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1.0  INTRODUCTION 


Dodge  (Reference  1**)  has  presented  a method  for  numerical  solutions 
to  the  momentum  and  continuity  equations  of  incompressible,  steady 
flow.  The  heart  of  this  method  is  the  substitution  of  Equation  (1) 
into  the  momentum  Equation  (2) : 


W = V<p  + u 


(1) 


DW 

Dt 


-VP/p  + V • ve(def  W) 


(2) 


where  s 


def  W = VW  + (VW) 


W = is  the  ensemble  averaged  velocity  vector 
P = the  ensemble  averaged  static  pressure 
p = the  density 

ve  ■ the  effective  kinematic  viscosity 


(VW) 


transpose  (VW) 


♦♦References  are  located  immediately  following  the  Table  of  Contents 
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The  result  is  Equation  (3)  for  u and,  with  the  application  of  continuity. 
Equation  (4)  for  potential. 

(W*V)  u + (u*  V)  V$  = veV  u + 7 (72tf>) 

+ (def  u)  Vve  + 2 ( V ( V4>  )]*  Vve 

V 2 = - V • u (4) 


Equation  (4)  is,  of  course,  elliptical  and  a direct  candidate  for 
relaxation  solutions.  During  the  relaxation  solution,  the  right-hand 
side  is  assumed  known  from  a previous  solution  to  Equation  (3) . Since 
Vp  is  directly  related  to  static  pressure,  Equation  (5)  , one  can  antici- 
pate that  its  variation,  unlike  W,  is  slow  near  the  wall.  Thus  the 
relaxation  grid  can  be  relatively  uniform  and  sparse. 

(V<{>*V)V<f>  = -VP/p 

The  result  is  the  relaxation  of  a single  dependent  variable  on  a 
sparse  grid  instead  of  four  or  more  dependent  variables  on  a dense 
grid.  The  disadvantage  is  the  necessity  of  solving  Equation  (3), 
and  iterating  between  Equations  (3)  and  (4). 

The  iteration  process  is  discussed  by  Dodge  (Reference  1).  The 
solution  to  Equation  (3)  would  be  easy  if  it  were  parabolic  in  a 
quasi-streamline  direction.  Under  these  conditions  it  could  be 
solved  by  a one-pass  marching  solution  analogous  to  those  commonly 
applied  to  boundary  layer  equations.  The  advantage  of  such  an 
approach  is  obvious,  since  a solution  is  very  rapid  even  when  a 
dense  grid  system  is  utilized.  However,  although  marching  solutions 
are  attractive,  several  obstacles  remain  for  their, direct  application 
to  Equation  (3).  Among  these  are  the  following: 
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o The  streamwise  diffusion  term  makes  Equation  (3) 

elliptic,  and  thus  creates  a non-well-posed  initial 
value  problem. 

o At  the  separation  point,  there  is  a possibility  of  some 
form  of  singular  behavior. 

o In  the  back  flow  region,  beyond  the  point  of  separation, 

the  flow  of  information  is  against  the  marching  direction, 
which  results  in  a non-well-posed  initial  value  problem  even 
without  streamwise  diffusion. 

The  following  sections  will  discuss  each  of  the  above  difficulties 
with  an  eye  toward  preserving  the  structure  of  the  numerical  method 
outlined  above.  The  final  section  presents  a comparison  between  a 
separated-f low  calculation  and  data  as  an  illustration  of  the  methods 
discussed  below. 
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2.0  STREAMWISE  DIFFUSION 


It  is  useful  to  examine  the  simpler  forms  of  Equation  (3)  to 
discover  pertinent  characteristics  of  this  equation.  Equation  (3), 
written  in  Cartesian  coordinates  x and  y with  constant  viscosity,  is 
given  by  Equation  (6) : 


3u 


3u. 


W 


+ W 


x 3x  "y  8y 

32u 


F U 


32u 


= v 


+ v 


3x 


3y 


32o>  . u lfjL_ 

x 3^2  Uy  8x3y 


+ v dr 


(6) 


The  general  similarity  between  Equation  (6)  and  the  boundary  layer 
equation  is  obvious.  If  the  boundary  layer  assumptions  are  applied 
to  Equation  (6),  Equation  (7)  results: 


W 


3u 

> 

3x 


W 


3u 

} 

y ay 


if* 

3x2 


32u. 


3y 


(7) 


This  transforms  directly  back  into  the  familiar  boundary  layer  equation 
by  the  substitution  of  Equation  (5)  for  the  potential  gradient.  It  is 
well  known  that  such  an  equation  is  parabolic;  however,  the  introduction 
of  streamwise  diffusion  changes  this  characteristic.  The  effect  of 
streamwise  diffusion  can  be  observed  by  considering  solutions  to  the 
following  two  linear  equations: 


3u 

7Z 


32u 

3y2 


0. 


(8) 


3u 

S3? 


32u 


3x 


3 u 


3y 


0. 


(9) 
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The  solution  to  Equation  (8)  is  given  by  Equation  (10): 


u 


eX2x/a 


+ 


e 


-Xyi 


(10) 


The  solution  to  Equation  (9)  is  given  by  Equation  (11): 


u = (e'ex  + Ae-ax)  (eXyi  + C2  e'^1)  (11) 

0 = J [l  + A + 4 xVc?] 
a = J [l  - A + 4 X2/o? 


For  a well-posed  initial  value  problem,  the  solution  at  x must  be 
bounded.  This  is  clearly  so  for  Equation  (8)  only  when  a is  negative. 
Otherwise,  there  exists  some  arbitrarily  large  frequency,  X,  that  will 
cause  an  arbitrarily  large  solution  at  x.  a being  negative  corresponds 
to  the  normal  viscous  flow  case  without  separation. 

Independent  of  a,  however,  Equation  (9)  cannot  be  represented  as 
an  initial  value  problem.  One  coefficient  or  the  other  (a  or  8)  is 
always  positive.  Dodge  and  Lieber  (Reference  2) , using  the  stability 
theory  of  Warming  and  Hyett  (Reference  3) , show  that  the  direct  conse- 
quence of  such  a non-well-posed  problem  is  that  no  consistent,  stable 
forward-marching  method  exists  using  finite-difference  operations.  In 
forward  flow,  the  number  of  steps  to  an  obvious  instability  strongly 
depends  on  the  magnitude  of  a.  However,  even  though  observable  insta- 
bilities may  not  be  present,  the  inclusion  of  streamwise  diffusion 
into  the  difference  operator  serves  no  purpose,  since  it  cannot  result 
in  increased  accuracy. 
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For  high  Reynolds  number  flows  where  streamwise  diffusion  is  small, 
the  obvious  assumption  is  to  neglect  its  effect.  This  assumption  has 
been  applied  beyond  the  boundary  layer  by  Dodge  (Reference  1)  and 
Pratap  et.  al.  (Reference  4) . In  the  boundary  layer  at  Reynolds  numbers 
such  that  the  boundary  layer  thickness  is  small  compared  to  axial 
dimensions,  neglect  of  the  effect  is  a very  good  assumption.  Figure  1 
shows,  as  an  example,  the  balance  of  a number  of  terms  in  the  momentum 
equation  for  a case  of  a linear  decrease  in  free-stream  velocity.  All 
terms  result  from  a boundary  layer  finite-difference  solution.  The 
streamwise  correction  was,  of  course,  not  used  in  the  differential 
equation,  but  calculated  as  an  after  result. 

Although  streamwise  diffusion  is  small  in  the  boundary  layer, 
neglect  of  it  may  not  be  justified  in  a separated  region.  A great 
deal  depends  on  the  length  scales  involved,  and  generalizations  are 
difficult.  If  the  length  of  a separated  bubble  is  similar  in  size  to 
its  width,  streamwise  diffusion  would  conceivably  play  a significant 
role.  Because  the  bubble  dimensions  are  strongly  influenced  by 
geometric  as  well  as  Reynolds  number  considerations,  a general  method 
for  solution  with  imbedded,  separated  bubbles  may  well  require  the 
inclusion  of  streamwise  diffusion  in  the  momentum  balance. 
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U = U -a  x 
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Figure  1.  Comparison  of  Terms  Through 
the  Boundary  Layer. 
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3.0  AT  THE  SEPARATION  POINT 

The  obvious  mathematical  similarity  between  the  marching  equation, 
Equation  (3),  and  the  boundary  layer  equation  suggest  that  perhaps 
a study  of  the  boundary  layer  equation  would  yield  insight  to  the 
performance  of  solutions  near  separation.  In  a recent  boundary  layer 
study.  Carter  and  Wornom  (Reference  5)  referred  to  a singularity 
existing  in  the  boundary  layer  solution  at  the  separation  point. 
Catherall,  Stewartson,  and  Williams  (Reference  6),  in  the  examination 
of  boundary  layers  with  normal  injection,  pointed  out  that  the  speci- 
fication of  a fixed  rate  of  injection  at  the  point  of  separation  leads 
to  an  inability  to  match  a prescribed  free-stream  velocity.  For  this 
case,  the  momentum  equation  at  the  wall  relates  first  and  second 
derivatives  of  the  velocity  profile  and,  by  inference,  all  other 
derivatives.  If  the  first  derivative  of  velocity  is  zero,  as  it  must 
be  at  a separation  point,  then  all  derivatives  are  zero.  Thus,  with- 
out some  form  of  singularity,  the  prescribed  free  stream  cannot  be  met. 
However,  with  no  injection,  this  form  of  singularity  does  not  exist. 

Catherall  and  Mangier  (Reference  7)  postulated  the  existence  of 
a singularity  at  separation  based  entirely  on  observed  numerical  insta- 
bilities. The  earlier  studies  of  numerical  methods  applied  around 
separation  were  carried  out  by  Leigh  (Reference  8) . Goldstein 
(Reference  9)  attempted  to  verify  the  existence  of  such  a basic  singu- 
larity by  asymptotic  expansions.  Stewartson  (Reference  10)  however, 
pointed  out  that  although  the  asymptotic  expansion  is  singular,  this 
provides  no  guarantee  that  the  equation  itself  is  singular.  Secondly, 
if  such  a singularity  exists,  it  is  so  narrow  as  to  preclude  detection 
by  normal  numerical  methods.  In  addition,  the  above  numerical  studies 
did  not  take  into  account  the  fact  that,  in  the  separated  bubble,  the 
boundary  layer  equation  is  a non-well-posed  initial  value  problem  and, 
without  modification,  should  not  be  soluble  by  any  f orward-marching 
method. 
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Still,  the  region  around  separation  is  extremely  critical  to  the 
numerical  technique  proposed  herein.  The  solution  proposed  for 
Equation  (3)  is  similar  to  that  of  the  boundary  layer  equation  in 
that,  through  <J>,  a prescribed  pressure  field  is  imposed;  however, 
it  is  different  in  that  continuity  is  solved  externally  to  the  marching 
equation.  Consequently,  a numerical  study  utilizing  the  boundary  layer 
equations  was  undertaken  to  observe  the  instabilities  in  and  around  the 
separation  point. 


The  solution  method  was  a basic  implicit  method  similar  to  that 
described  by  Schlichting  (Reference  11) . The  free-stream  velocity 
was  prescribed  as  linearly  decreasing  with  distance  along  the  plate, 
x.  Ahead  of  separation,  the  solution  proceeds  in  a normal  fashion. 
Downstream  of  separation,  however,  the  following  assumptions  were 
introduced  in  the  back  flow  region  in  an  attempt  to  produce  a stable 
solution. 


Case 

A. 

No 

change 

Case 

B. 

No 

convection 

parallel 

Case 

C. 

No 

convection 

parallel 

Case 

D. 

No 

convection 

parallel 

set  to  zero 

Case 

E. 

Normal  velocity  set  to 

to  the  wall 

or  normal  to  the  wall 

to  the  wall.  Normal  velocity 

zero. 


The  above  assumptions  only  apply  to  the  regions  where  backflow 
occurs.  Instabilities  were  not  noted  until  they  became  obvious.  As 
discussed  in  the  previous  section.  Case  A should  be  unstable.  The 
elimination  of  free  stream  convection,  B,  removes  the  obvious  problem 
of  the  momentum  equation  being  non-well-posed.  However,  the  normal 
component  of  velocity  still  has  upstream  influence  through  the  con- 
tinuity equation.  The  elimination  of  both  normal  and  parallel  con- 
vection should  completely  eliminate  the  difficulty  with  the  momentum 
equation.  The  results  for  the  first  three  cases  (A,  B,  and  C)  were  an 
obvious  instability  by  about  10  percent  (of  the  length  to  the  separation 
point)  downstream  of  the  separation  point. 
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Removing  the  objectionable  terms  from  the  momentum  equation  had 
little  noticeable  effect  on  the  results.  Case  B,  when  applied  to 
Equation  (3)  by  Dodge  (Reference  1)  did,  however,  eliminate  all 
observed  numerical  instabilities.  The  obvious  difference  is  that  the 
boundary  layer  equations  include  a continuity  equation  with  upstream 
influence  where  the  equations  of  Dodge  do  not.  Consequently,  Cases  D 
and  E were  tried.  Case  D eliminates  the  objectionable  continuity  and 
parallel  convection  terms.  Case  E is  included  as  a check  case  in 
which  continuity  is  removed  but  parallel  convection  is  not.  The 
results  were  as  follows: 

Case  D:  No  obvious  instability  within  at  least  100  percent  of 

the  separation  length  downstream. 

Case  E:  Obvious  instability  at  about  10  percent  downstream 

of  the  separation  point. 

The  success  of  Case  D indicates  that  if  the  numerical  problems 
created  by  backflow  can  be  overcome,  a solution  for  an  arbitrarily 
imposed  potential  and  thereby  pressure  gradient  is  obtainable.  Ignor- 
ing parallel  convection  is,  however,  a suitable  assumption  only  if  the 
backflow  region  is  sufficiently  small.  Streamwise  convection  can 
rapidly  become  significant  as  indicated  by  Figure  2.  Streamwise  con- 
vection is  calculated  from  the  results  of  Case  D.  Even  though  the 
point  in  question  is  only  slightly  downstream  of  the  separation  point, 
free-stream  convection  can  hardly  be  ignored.  This,  then,  leads  to 
the  numerical  methods  described  in  the  following  section. 
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Figure  2.  Comparison  of  Terms  by  Momentum  Balance 
Through  the  Boundary  Layer  Method  of 
Solution  D. 
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4.0  THE  BACKFLOW  REGION 

As  already  observed  in  Section  II,  the  effect  of  backflow  is  to 
assure  that  Equation  (3)  cannot  be  represented  as  a well-posed  initial 
value  problem.  Returning  to  the  linear  Equation  (9) , a technique  can 
be  developed  to  circumvent  this  problem.  The  dependent  variables  are 
replaced  by  Equation  (12) : 


u(x,y)  = C (x) e lwy 


(12) 


This  can  be  done  without  loss  of  generality  since  all  frequencies  can 
be  represented  by  Equation  (12) . The  result,  for  constant  a,  is  an 
ordinary  differential  equation  for  the  coefficient,  C(x): 


- w2C  =-  0 . 


(13) 


The  solution  to  Equation  (13)  is  given  by  Equation  (14): 


, , . -ax 

C (x)  = A exp  l — 2” 


|l  + / 1 + 4u)2/a2]  ) 


+ B exp  ( — ^ | 1 “ / 1 + 4ui^/a"2  | ) 


(14) 


Note  that  this  equation  is  similar  to  Equation  (11) . However,  Equa- 
tion (14)  does  not  violate  the  conditions  of  a well-posed  initial 
value  problem,  because  its  solution  depends  in  a continuous  and  bounded 
fashion  on  the  initial  conditions  (a)  and  (b) . Following  the  nomencla- 
ture of  Richtmyer  and  Morton  (Reference  12) , if  the  solution  to  an  ini- 
tial value  problem  can  be  represented  by  a transformation  given  by 
Equation  (15) , then  the  following  conditions  are  required  for  the  prob- 
lem to  be  well  posed: 
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(a) 

(b) 


Eq(x)  is  continuous 

Eq(x)  is  uniformly  bounded — i.e., 

||  E (x)  ||  < K for  x . < x 

° min  — 

where  K must  be  a finite  constant 


< 


xmax 


u(x) 


UoEo(x) 


(15) 


Equation  (9)  violates  condition  (b) , whereas  Equation  (13)  on  the  other 
hand  does  not  violate  condition  (b) . Only  a single  finite  frequency  is 
being  considered;  thus,  Eq(x)  is  uniformly  bounded,  and  stable  forward- 
marching, finite-difference  operations  consistent  with  Equation  (13) 
can  be  found. 


This  leads  directly  to  consideration  of  such  a technique  applied 
to  the  more  complicated  non-linear  equations  of  marching.  Equation  (3) 
expressed  in  a two-dimensional  curvilinear  coordinate  system  can  be 
represented  by  Equation  (16) : 


A u + 
o 


A ~ 

nj^  «'v 


— + A 3u 

axx  + *2  B5J 


+ A 


11 


32u 
3x,  : 


+ A 


~ 2 
9 u 

22  a 2 
0X2 


s = 0 


(16) 


In  this  case,  only  one  of  the  dependent  variables  is  shown,  u represent- 
ing either  ux  or  u2,  depending  on  the  situation.  Additionally,  the  co- 
efficients of  the  equation  are  non-linear  functions  of  u^,  u2,  x^j^ , and 


The  dependent  variable  can  be  represented  by  an  infinite  series  of 
a sufficiently  general  set  of  functions: 


“‘V**1  - c0(x1)g0(x  > + Z cn(x  ) gn(x  ) 

n=  1 


(17) 
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Substitution  of  Equation  (17)  into  (16)  yields  Equation  (18): 

E A C q + A,C'  g + A-C  q'  + A C'a  + A C a"  I 

n=l  \ ° n n 1 n yn  2 nyn  Mll'"ngn  A22  ngn  / 


-S-A  C a 
o O^O 


-aic;  g0 


- A2cog; 


- AnC"g 
1 1 o^o 


A~~C  g 
2 2 o^o 


(is: 


Treating  the  zeroth  function  separately  allows  it  to  be  chosen  to  match 
exactly  the  boundary  conditions  resulting  in  simpler  numerical  solu- 
tions . 

During  the  course  of  the  study,  several  methods  for  obtaining  the 
unknown  coefficients  were  attempted.  Dodge  and  Lieber  (Reference  ?)  de- 
scribe these  developments.  Two  are  noteworthy  enough  to  be  described 
herein. 


The  first,  which  will  be  referred  to  subsequently  as  the  least 
squares  method,  represents  only  an  approximate  solution  to  the  unknown 
coefficients.  Each  set  of  terms  in  Equation  18  may  be  equated  to  a 
function  en(x1/  x2)  representing  an  error,  as  follows: 


A C + A.C  ' 
on  ln 


+ A?cn9'/‘Jn  + AnC^  + Aoo  c S"/g 

<Jnn  n 11  n 22  nJn^n 


+ d 

n 


e 

n 


<x! ,X2  > 


(19) 


Then  the  term  is  given  by: 


s 


E 

n=l 


d g 
n^n 


d g 
o^o 


If  the  right-hand  side  was  identically  zero  for  all  n,  then  an  exact 
solution  would  be  obtained.  Note,  however,  that  this  is  not  the  only 
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way  a solution  can  be  reached.  The  other  way  is  for  e to  be  non-zero 

n 

for  a particular  value  of  n,  but  canceled  by  other  frequency  components. 
The  type  of  solution  depends  on  the  series  choice  and  the  problem  type. 

If  the  error  en  is  approximated  at  a discrete  number  of  points  M, 
in  the  X2  direction,  then  a total  squared  deviation  parameter  is  de- 
fined by  Equation  (20) : 


E 

n 


M 

E 

m=l 


(20) 


The  solution  nearest  the  identically  zero  solution  could  be  ob- 
tained by  minimizing,  En,  the  sum  of  the  squares.  To  do  this,  the  un- 
known coefficients  are  expanded  in  difference  equations  as  indicated 
by  Equation  (21): 


lll 


Ax, 


c (k)  _ 2C  (k-l)  (k-2) 

n n n 


"1 

2Ax, 


3 C (k)  - 4 C (k_1)  + r (k‘2) 
n n n 


where : 


(k) 


n 


A_f' 
2 n 


/f  + A 


22 


f" 

n 


/f  + A, 


n 


+ d = e 


n 


(21) 


xl  = xl  + kAx  ^ 
o 

Collecting  terms  in  Equation  (21),  a simple  relation  for  the  unknown 

coefficient  C ^ is  obtained, 
n 
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B C (k)  - P = e 

m n m n 


where : 


B = 
m 


A,  , /Ax^ 2 + ^1  + A g'/g  + A0.g"/g 

x 2Ax  1 n n 


'll' 


+ A 


(22) 


P 


m 


c <k-x) 

n 


2A 


11  + 


c (k-2> 


+ 


d (k> 


Differentiation  of  Equation  (20)  by  the  unknown  coefficient,  C , 

(k)  n 

results  in  an  evaluation  of  C yielding  the  least  square  error: 


M 

Z 

m=l 


P_B 
m m 


M 

Z 

m=l 


(23) 


Solutions  to  Equation  (24)  were  obtained  utilizing  this  technique: 


(W  +u) 

00 


9u 

9x 


v 


subject  to  the  boundary  conditions: 

u (x , 0)  = -w 

00 


2 ‘ 

3 u , ,, 

— r (x,6)  = 0. 

3x^ 


(24) 


Since  this  is  the  equation  for  flow  over  a flat  plate  with  no  pressure 
gradient,  backflow  solutions  were  obtained  by  marching  downstream  and 
then  reversing  the  process  in  an  attempt  to  recover  initial  conditions. 
Simultaneous  solutions  were  made  with  conventional  finite-difference 
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methods.  The  finite-difference  solution  is  compared  to  the  least 
squares  solution  for  the  downstream  marching  in  Figure  3.  In  the  re- 
verse direction,  finite-difference  solutions  become  obviously  unstable 
within  a few  marching  steps.  Figure  4 displays  the  results  of  the 
least  squares  calculation.  Note  that  the  return  solution  is  a close 
approximation  to  the  initial  conditions.  Figure  5 shows  the  variation 
of  spectral  coefficients  involved  in  this  process.  The  series  used  is 
given  by  Equation  (25)  : 


N 

u(x,y)  = -W^  + E cn  (x)  sin  ((2n-l)y  tt/6 ) ) (25) 


Note  that  the  boundary  conditions  are  automatically  preserved  by  such 
a choice.  Secondly,  it  is  useful  to  observe  the  low  number  of  func- 
tions actually  utilized.  It  might  be  expected  that  such  a limited  set 
could  not  describe  a turbulent  boundary  layer.  However,  the  number  of 
functions  depends  on  which  functions  are  used,  and  the  orthogonal 
transformation  applied  to  the  grid  system.  As  will  be  seen  in  a later 
section,  good  results  can  be  obtained  with  only  a modest  number  of 
coefficients,  even  for  a very  complicated  situation. 

Such  a least  squares  method  is  suitable  for  obtaining  insight 
into  the  method,  and  showing  its  capability  to  solve  backflow  problems. 
However,  it  is  not  practical  for  more  complicated  situations  involving 
pressure  gradients  and  curvilinear  geometries.  Dodge  and  Lieber 
(Reference  2)  show  that  the  simple  addition  of  a pressure  gradient  to 
Equation  (24)  results  in  large  solution  errors  by  the  least  squares 
method.  Consequently,  the  following  method  was  developed. 
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Figure  5.  Spectral  Coefficients. 
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5.0  METHOD  OF  WEIGHTED  RESIDUALS 


A better-based  method  of  solution  can  be  obtained  by  multiplying 
Equation  (18)  by  one  member  of  the  functional  set,  g^,  and  then  inte- 
grating over  the  solution  region. 


n=l  jC"  f6(A°9ngk  + dx2 


+ cn  'j  A1  Vkdx2  + Cnf 


A. , g g,  dx_ 
11  ^r.^k  2 


6 /*i5  g 

= - [ s^kdx2  - CoJ  A1  ^o^kdx2-Co  / AllVkdX2 

Jo  o Jo 

~Co[(A22  gogk  + Aogogk  + A2go  gkfX2 


(26) 


The  result,  for  a finite  number  of  functions  is  a set  of  coupled  ordi- 
nary differential  equations  for  the  unknown  coefficients,  Cn<  If  g^ 
and  its  derivatives  are  orthogonal  over  the  interval  (0,6),  and  if  the 
coefficients  of  the  differential  equation  are  not  a function  of  x2, 
the  off-diagonal  terms  become  zero  and  the  equations  decouple.  The 
linear  equation  solution  given  by  Equation  (14)  then  falls  out  immedi- 
ately. In  general,  with  a non-linear  curvilinear  transformation  to 
Equation  (3) , the  off-diagonal  terms  of  Equation  (26)  are  not  large  and 
a reasonably  non-singular  matrix  results.  This  technique  is  similar  to 
several  methods  reported  by  Ames  (Reference  13)  for  fluid  dynamic  prob- 
lems . 


The  choice  of  function  set  can  be  somewhat  artistic.  The  func- 
tional set  was  chosen  so  that  the  anticipated  boundary  conditions  were 
met  automatically.  These  conditions  on  the  wall  are  given  by  Equations 
(27)  and  (28): 


a 
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u2  (xx,0)  = 


1 94>  (xlfo) 

h„  3x, 


-V  (x2) 


(27) 


u-^x^o)  = 1/h^ 


dtp  ( X^O ) 
3xi 


(28) 


Either  a zero  value  of  u or  a zero  derivative  of  u was  applied  to 
the  outer  boundary  depending  on  the  series  used. 

The  following  functional  sets  were  examined. 


Set  A - Sine  Series 


W 


. I 


) sin  X x- 
l n 2 


n=0 

n^O 


X = ( 2n-l)  n/6 
n 


Set  B - Modified  Chebyshev  Series 


g (x2) 


Tn  + lml  dk  Tn+k 


where: 


2x2  -1.0 

~T 


Tn(y*)  = Chebyshev  polynomial  of  order  n 


the  coefficient  d^  is  chosen  so  that 
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Set  C 


gQ(o)  = 1.0 
gn(o)  =0  n _> 

gn(6)  = 0 n >_ 

g ' (6)  = 0 n 1 

Sine-Cosine  Series 


0 

0 

0 


gn  (x2)  = sin  Anx2 

\ nir 

n = T 


n > o 


go  (x2)  = COS  77X2 


Set  D - Poly  Sine  Series 


W = (1_x2/6) 


gn(x2>  = sin  (nvx^/S)  n > o 

Each  of  the  above  methods  were  compared  against  solutions  to  Equation 
(29): 

9x  dx2  9y  129) 

where: 

" - & * ■ 

The  boundary  layer  type  equation  corresponding  to  Equation  (29)  is 
given  by  Equation  (30) : 

. 1301 
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A Pohlhausen- type  solution  was  used  for  comparison.  The  Pohlhausen 
solution  differs  from  that  of  Schlichtiny  (Reference  11) , however,  in 
that  no  continuity  equation  is  present.  The  development  of  this  par- 
ticular Pohlhausen  solution  is  given  by  Dodge  and  Lieber  (Reference  2) 
and  will  not  be  included  here.  Various  amounts  of  diffusion  and  ac- 
celeration were  applied.  Initial  conditions  were  based  upon  a 
Pohlhausen  profile  at  a very  small  value  of  x. 

Since  the  Pohlhausen  solution  is  not  exact,  it  cannot  be  compared 
directly,  particularly  for  diffusing  cases  where  the  Pohlhausen  solu- 
tion predicts  a larger  displacement  thickness  than  actually  exists. 

The  Pohlhausen  solution  is  then  used  as  verification  of  general  trends 
and  internal  comparisons  must  be  made  to  differentiate  somewhat 
imperfectly  between  methods. 

The  modified  Chebyshev  series  illustrates  the  importance  of  the 
orthogonality  or  near-orthogonality  in  the  function  set.  The  motiva- 
tion behind  this  set  of  functions  rests  in  the  fact  that  Chebyshev 
polynomials  can  be  evaluated  much  more  rapidly,  due  to  the  existence 
of  an  algebraic  recursion  relationship  between  successive  Chebyshev 
polynomials.  However,  as  a plot  of  the  first  three  functions  shows 
(Figure  6) , the  functions  are  far  from  orthogonal.  They  are  suffi- 
ciently non-orthogonal  that  the  resulting  matrix  is  dif f icult-to- 
impossible  to  invert  and  no  solutions  were  obtained. 

The  original  choice  of  the  function  set  was  Set  A.  Set  C was 
chosen  since  the  cosine  functions  more  closely  approximate  a boundary 
layer  profile  (varying  from  unity  at  *2=0  to  zero  at  x2=l^  than  does 
the  constant  of  Set  A.  Figures  7,  8,  and  9 show  the  results  for  ac- 
celerating flow,  and  two  points  in  decelerating  flow.  Figure  9 con- 
tains several  solutions  with  differing  spacing  in  an  attempt  to  ex- 
plore the  convergence  of  the  method.  From  the  baseline,  the  step  size 
in  the  x direction  and  then  in  the  y direction  was  reduced.  Higher  ac 
curacy  solutions  should  proceed  from  squares  to  diamonds  to  triangles. 
It  is  not  unreasonable  to  see  such  a trend  toward  higher  accuracy,  as 
indicated  in  Figure  9. 
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Figure  7.  Comparison  of  Calculated  Velocity  Profiles  for 
Accelerating  Flow. 
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Figure  9.  Comparison  of  Calculated  Velocity  Profiles. 
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The  poly-sine  series  (Set  D)  was  chosen  to  correct  one  deficiency 
in  the  sine-cosine  series  (Set  C) . The  cosine  has  the  wrong  inflec- 
tion across  the  boundary  layer.  The  polynomial  of  Set  D does  not. 
Figures  10  and  11  show  the  results  for  accelerating  and  decelerating 
flow  compared  to  Pohlhausen  and  Set  C.  Little  significant  difference 
between  individual  solutions  was  observed.  However,  the  magnitude  of 
first-  and  second-order  coefficients  sharply  reduced  between  Sets  C 
and  D. 

The  sine  series  (Set  A)  was  the  original  series  (on  a chronologi- 
cal basis) . Figures  12  and  13  show  comparisons  between  calculation 
methods  for  accelerating  and  decelerating  flow.  Again,  little  differ- 
ence is  observed.  However,  higher  order  coefficients  are  largest  in 
magnitude  with  Set  A. 

The  conclusions  of  the  above  study  were:  first,  few  differences 

between  orthogonal  function  based  methods  (Set  A,  C,  and  D)  were  ob- 
served; second,  it  is  important  to  have  function  sets  that  are  orthog- 
onal, or  nearly  so;  and  third,  the  general  solution  agrees  with  expec- 
tations. With  this  experience,  a program  to  solve  for  the  flow  about 
an  isolated  body  in  a wind  tunnel  without  sidewall  boundary  layers  was 
undertaken. 
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Figure  10.  Comparison  of  Calculations  for 
Accelerating  Case. 
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Figure  11.  Comparison  of  Calculations 
for  Decelerating  Case. 
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Figure  12.  Comparison  of  Calculations  for 
Accelerating  Case. 
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6.0  COMPARISON  TO  DATA 

A demonstration  of  the  above  method  was  undertaken  for  a full- 
scale  problem.  Comparison  was  made  to  the  experimental  data  taken  by 
Schubauer  and  Klebanoff  (Reference  14) . The  configuration,  consisting 
of  a very  long,  large  shaped  body  in  a wind  tunnel,  is  shown  in  Figure 
14,  along  with  some  predictions. 

The  series  used  over  the  body  was  Series  A.  A complete  Fourier 
series  was  used  ahead  of  and  behind  the  body.  Transition  to  the  body 
and  to  the  wake  requires  some  attention.  In  a difference  equation  ap- 
proach to  the  solution,  the  nodes  near  the  leading  and  trailing  edges 
represent  only  a local  discontinuity.  The  coefficients  of  the  series, 
however,  represent  global  profile  information  and  are  discontinuous  at 
the  leading  edge.  Coefficients  at  this  point  must  be  obtained  by  inte- 
grating the  upstream  profile  to  a discontinuous  profile,  one  that  goes 
to  no-slip  at  the  location  of  the  wall. 

The  grid  system  utilized  results  from  a numerical  transformation. 
Such  a grid  system,  without  the  closely  spaced  points  near  the  wall,  is 
shown  in  Figure  15.  Wall  points  are  placed  at  the  users  discretion. 
Good  results  have  been  obtained  using  a logarithmic  spacing  with  a fac- 
tor of  several  thousand  between  the  smallest  and  largest  increments. 
Something  on  the  order  of  20  points  are  used  across  the  boundary  layer 
with  a like  number  across  the  rest  of  the  tunnel  passage.  The  relaxa- 
tion solution  was  carried  out  on  a sparse  subset  of  the  marching  grid 
utilizing  its  own  transformation  coefficients. 

The  comparison  between  experiment  and  data  in  incompressible  flow 
is  sharply  complicated  by  the  need  for  a turbulence  model.  Needless 
to  say,  no  agreement  on  an  adequate  model  exists,  especially  when  com- 
plicated by  separated  bubbles.  As  such,  the  comparison  can  only  be 
used  as  further  evidence  of  the  reasonableness  of  the  results.  To 
this  end,  the  basic  turbulence  model  was  essentially  that  utilized  by 
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Figure  14.  Shaped  Body  in  Wind  Tunnel. 
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CIRCULAR  CYLINDER  IN  CROSS  FLOW  RA0IU8-1.8 


Figure  15.  Full  Relaxation  Grid. 
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Dodge  (Reference  1)  successfully  in  three-dimensional  flows.  This 
basic  model  is  a mixing- length  model  with  the  mixing  length  modified  by 
van  Driest' s (Reference  15)  correction  for  wall  shear  and  the  Cebeci- 
Smith  (Reference  16)  correction  for  pressure-gradient  effects  on  tur- 
bulence. The  mixing  length  increases  until  a cutoff  value  is  reached. 
This  cutoff  is  a function  of  free  stream  turbulence.  Dodge  reports 
a cutoff  at  approximately  0.02  of  chord  for  levels  of  turbulence  con- 
sistent with  internal  flows.  Comparisons  to  data  were  made  at  the 
trailing  edge,  and  undoubtedly,  significant  error  occurred  near  the 
leading  edge.  Although  these  values  of  £max  are  somewhat  larger  than 
those  suggested  by  Launder  and  Spalding  (Reference  17)  they  are  repre- 
sentative of  the  values  required  to  predict  a fairly  wide  range  of 
internal  flows.  The  three  following  approximations  were  utilized  for 
calculations  on  the  Schubauer  and  Klebanoff  configuration. 

(1)  A constant  cutoff  of  £ = 0.585  feet 

max 

(2)  A variable  cutoff  with  two  different  distributions 
(see  Figure  16) . 

(3)  An  Alber  cutoff  on  the  value  of  kinematic  viscosity 
at  the  edge  of  the  boundary  layer. 

u = 0.0168  v <5* 
e e 

where:  Wg  is  the  edge  velocity 

6*  is  the  displacement  thickness. 

Results  for  all  three  models  are  shown  in  Figure  17  for  displace- 
ment thickness  and  Figure  18  for  shape  factor.  As  expected,  the 
constant  cutoff  produced  a too-rapid  early  growth  with  a too-slow 
growth  near  the  trailing  edge. 
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Figure  17.  Comparison  of  Experimental  and  Calculated 
Displacement  Thicknesses. 


74-211196(6) 
Page  39 


AIRESEARCH  M AN  U FACTU  Rl  N G COMPANY  OF  ARIZONA 

A Q>V««<0*  O*  CO*#D*  >OM 

PHOENIX.  ARIZONA 


Figure  18.  Comparison  of  Experimental  and  Calculated 
Shape  Factor. 
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The  variable  cutoff  does  quite  well,  however.  The  Alber  model  suffers 
in  a full  Navier-Stokes  solution  with  the  need  to  define  such  nebulous 
quantities  as  the  freestream  velocity  and  the  displacement  thickness. 
The  numerical  method  used  resulted  in  an  apparently  spurious  inter- 
action between  cutoff  values  and  solutions,  yielding  the  poor  results 
shown.  Undoubtedly,  it  was  a problem  that  could  be  overcome,  but  it 
was  judged  not  to  be  worthwhile  to  pursue  it  further.  Figure  19  shows 
the  separation  points  of  the  various  techniques  and  a calculated 
separated  bubble.  An  earlier  separation,  closely  confined  to  the 
surface,  is  predicted  by  Model  (2) ; however,  the  bubble  for  Model  (2) 
has  a reasonable  position  and  size  when  compared  to  data. 


Velocity  profiles  are  shown  for  the  following  cases  to  illustrate 
the  effect  of  successive  passes  through  the  relaxation  equation  (Equa- 
tion 14)  with  various  flavors  of  the  first  two  turbulence  models. 


Case  A 1 outer  loop  (one  relaxation,  one  march) , 

Model  1 with  l =0.585 
max 

Case  B 1 outer  loop.  Model  2 with  ^max  variable  as 
specified  by  Figure  16. 

Case  C 1 outer  loop.  Model  2 with  l variable  as 

max 

specified  by  Figure  16. 

Case  D 3 outer  loops  (3  relaxations,  3 marches), 

Model  2 with  ^max  variable  as  specified  by 
Figure  16. 

Only  the  calculation  points  in  the  outer  boundary  layer  portion  of  the 
flow  are  shown.  Additional  points  are  closer  to  the  wall  and  out  in 
the  channel.  Figures  20,  21,  22,  23,  and  24  compare  profiles  at  dis- 
crete locations  along  the  body.  As  can  be  expected.  Case  D offers  a 
reasonable  fit  to  data.  The  effect  of  iteration  is  most  pronounced 
near  the  trailing  edge  where  the  separation  region  changes  from  itera- 
tion to  iteration,  indicating  the  importance  of  the  static  pressure  in 
this  region. 
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Figure  19.  Experimental  and  Calculated  Separation  Points. 
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Figure  20.  Comparisons  of  Calculations  to  Experimental  Boundary 
Layer  Velocity  Profile  at  X = 3.5  Ft. 
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Figure  21.  Comparisons  of  Calculations  to  Experimental  Boundary 
Layer  Velocity  Profile  at  X = 10.5  Ft. 
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Figure  22. 


Comparisons  of  Calculations  to  Experimental  Boundary 
Layer  Velocity  Profile  at  X = 17.5  Ft. 
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Figure  23.  Comparisons  of  Calculations  to  Experimental  Boundary 
Layer  Velocity  Profile  at  X = 21.0  Ft. 
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Figure  24.  Comparisons  of  Calculations  to  Experimental  Boundary 
Layer  Velocity  Profile  at  X = 25.77  Ft. 
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The  results  of  these  calculations  were  not  intended  to  provide  a 
turbulence  model  of  any  great  universality,  but  rather  to  demonstrate 
the  viability  of  this  technique  for  solution  of  the  fluid  dynamic  equa- 
tions. Indeed,  the  effect  of  turbulence  complicates  the  evaluation 
since  the  base  accuracy  of  the  method  cannot  be  assessed.  Neverthe- 
less, the  comparison  is  good  with  a model  that  has  some  reasonable 
foundation. 
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7.0  CONCLUSIONS 


The  above  sections  have  outlined  a method  that  offers  the  poten- 
tial for  rapid  solutions  to  the  Navier-Stokes  equations.  Run  times 
for  the  problem  described  in  Section  6.0  were  less  than  7 minutes  on  a 
CDC  6400.  Substantial  improvement  to  this  could  easily  have  been  made 
by  improving  the  coding  efficiency  in  both  the  relaxation  and  marching 
sections.  The  solution  method  reduces  the  expensive  relaxation  portion 
of  the  calculation  to  just  that  needed  to  compute  static  pressure. 

j t — k t $ 

The  vir&eotis  equation  is  solved  by  a single-pass  marching  process, 
which  overcomes  the  mathematical  difficulties  imposed,  by  transforming 
the  solution  in  the  direction  normal  to  the  flow  into  a spectral  plane. 
The  method  has  been  verified,  as  far  as  possible.  However,  a precise 
determination  of  accuracy  is  clouded  by  the  uncertainties  introduced  by 
turbulence  and  the  modeling  thereof.  Nevertheless,  it  has  been  shown 
that  the  method  yields  reasonable  results  both  in  a full  solution  to 
the  flow  around  a shaped  body,  and  in  a variety  of  plate  flows  with 
and  without  pressure  gradients,  Further  exercise  of  the  method  would 
be  greatly  enhanced  by  the  inclusion  of  compressibility  effects,  which 
would  allow  study  of  shock  boundary  layer  interactions  where  signifi- 
cant laminar  data  exists. 
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